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(57) Abstract 



Molecular modeling is performed using atomic parameters 
which include an anisotropic dipole polarizability tensor. Per- 
manent atomic multipole parameters may also be included in the 
model. Energy evaluations including contributions from polariza- 
tion energy and multipole interactions may be conducted which 
are useful in characterizing molecular properties for dmg discov- 
ery, materials evaluation, and other applications. 
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MODELING INTERACTIONS WITH ATOMIC PARAMETERS INCLUDING ANISOTROPIC DIPOLE POLARIZABILITY 



35 



_. , , , ^ . Packnrotind 

Field nf the ttwftitjnn- 

The invention relates to computation of interatomic and intermolecular interactions. Embodiments of the invention 
are applicable to chemical. physicaL and biological research, and to the development of new phamiaceutical compounds. 
Descrintinn nf ^|,e Ralatprt flff- 

Over the last few years, significant advances have been made in predicting, studying, and quantifying the nature 
of interatomic and intemiolecular interactions with'the use of computer simulations. Although from a purely scientific point 
of view, such computer simulations can be useful in testing and validating a theory conceming the nature of these 
interactions, computer simulations find especially useful application in reducing the time required to develop new materials 
with desirable properties. Materials which may be developed with additional efficiency through the use of computer 
simulations include polymers, pesticides, herbicides, pharmaceuticals, semiconductor materials for integrated circuits, and 
petrochemicals to name a few. 

Several modeling techniques are used in these environments. Typically, the model selected provides a user of the 
model with a particular compromise between physical accuracy and the computing resources required to run the model. 
initio quamum mechanical calculations can be perf omied with a high degree of accuracy, but are very expensive in terms of 
the computer time required to perform them. In those fields described above, where new polymers, drugs, etc. are being 
developed, it is more useful to use modeling techniques which require less imrestment in computing resources, so that more 
candidate materials or molecules can be analyzed in a shorter time period for the properties desired. 

Thus, it has, become common to model interactions between groups of atoms, molecules, proteins, and other 
structures by defining an atom to atom potential which acts between the atoms of the system being analyzed. Generally 
speaking, the atom-atom potential defines the energy of the atomic system as that energy varies with the coordinates of the 
atoms. Intramolecular "bonded" interactions may include temis defining energy as a function of bond lengths, bond angles, 
torsion angles, and out-of plane coordinates. Intemiolecular, or non-bonded potentials, typically include van der Waals 
interactions and electrostatic interactions. The benefit of these force field models is that the behavior of the atoms in the 
model is calculated using classical mechanics and electrostatics, which is significantly simpler computationally than 
performing the more mathematically complex quantum mechanical calculations. 

However, because the systems being analyzed do not in fact behave classically, the models include parameters 
associated with the force field terms which are selected to fit known quantum mechanical molecular and atomic interactive 
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behavior. In this way, a classically formulated force field is used to approximate quantum mechanical behavior at the 
atomic and molecular scale. A variety of force field models are known. Force field models are provided, for example, in U.S. 
Patent No. 5,612,894 to Wertz, and in -OREIDING: A Generic Force Field for Molecular Simulations.' J. Phys. Chm. 94. 
8897-8909 (1990). An additional force field model nicknamed TFF', and which was developed by several of the inventors 
of the subject matter of the present application, is described in^. Comp. Chem 15. 162-182 11994). and in ^. Am. Chem. 
Soc. 116. 2515-2525 (1994). The disclosure of U.S. Patent No. 5,612.894 and the Journal articles described above are 
hereby incorporated by reference in their entireties. 

It can be appreciated that the selection of appropriate parameters and force field functional dependencies is a 
significant factor in the success of the model. Furthermore, the number of fitted parameters used in the model relative to 
the number of measurable or ab initio calculable values relevant to the system being modeled should be as low as possible. 
A model with as many fitted parameters as observables has little predictWe value for systems which were not used in 
initially creating the fitted parameters. 

Currently, most force fields treat interatomic electrostatic interactions using a partial charge model in which each 
atom is assigned a net charge and Coulomb's Law is used to calculate forces on each atom due to the other atoms of the 
system. The DRIEDING and CFF force fields mentioned above are examples. Another alternative which has been devised is 
to model a molecule as a set of bond centered dipoles. Neither treatment adequately models the interaction between atoms 
and the local electric fields. Accordingly, new force field parameterization schemes are needed to increase agreement with 
experiment, to maximize the number of observables relative to the number of parameters, and to limit the necessity of 
performing computationally expensive calculations. 

i?limm'"^V 'he Invention 

The invention includes methods of evaluating a candidate molecule for suitability for a particular purpose. In one 
embodiment, the method includes selecting a candidate molecule and calculating a dipole moment induced in a first atom of 
the candidate molecule from a local electric field, wherein the induced dipole moment may be non-parallel to the local electric 
field. Using the calculated dipole momem. one or more physical properties of the candidate molecule may be predicted. 

In another embodiment, a method according to the invention comprises parameterizing electrostatic behavior of the 
candidate molecule with a plurality of atomic parameters associated with at least one atom of the molecule, wherein the 
plurality of atomic parameters includes elements of an anisotropic atomic dipole polarizability tensor. The method also 
includes determining a dipole moment induced in the atom due to a local electric field using the atomic dipole polarizability 
tensor and predicting one or more physical properties of the candidate molecule using the Induced dipole moment. 

Embodiments of the invention also include apparatus for modeling the geometry and energy of interaction between 
first and second groups of atoms. The apparatus may comprise a memory storing an anisotropic atomic dipole polarizability 
tensor for at least one of the first group of atoms. Also provided may be a processor for (1) modeling an electric field 
produced at least in part by the first and second groups of atoms. (2) retrieving the anisotropic dipole polarizability tensor, 
and (3) calculating a dipole moment induced in an atom of the first group by the electric field, and (4) calculating an 
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interaction energy between the first and second groups of atoms which includes a contribution from the induced dipole 
moment. The apparatus may further include an output device for reporting the calculated interaction energy. 

Another embodiment of the invention is a computer readable media having stored thereon commands which cause 
a general purpose computer to perform a method of modeling interactions between a first group of one or more atoms and a 
second group of one or more atoms. In one embodiment, the method comprises retrieving, for at least one of the first group 
of atoms, elements of an anisotropic atomic dipole polarizability tensor from a data storage device, modeling an electric field 
which is dependent on a relative orientation and separation of the first group of atoms and the second group of atoms. The 
method further includes calculating a dipole moment induced in the atom of the first group from the electric field using the 
anisotropic polarizability tensor, calculating an interaction energy between the first group of atoms and the second group of 
atoms which includes a contribution from the induced dipole moment, and outputting the interaction energy. 

Brief Description of the Drawings 
FIG. 1 is a perspective view of interacting atoms in a interaction model. 

FIG. 2 is a flowchart illustrating the development and use of atomic parameters for modeling electrostatic 
interactions according to one embodiment of the invention. 

FIG. 3 is an illustration of a water molecule orientation which may be used in calculating a portion of an atomic 
parameter set. 

FIG. 4 is a flowchart of one method of implementing atomic parameterization according to the invention. 

FIG. 5 is a flowchart of a geometry optimization method utilizing an embodiment of the invention. 

FIG. 6 is a flowchart of one method of calculating polarization energies in the geometry optimization method of 

20 Figure 5. 

FIG. 7 is a flowchart of one method of calculating polarization forces in the geometry optimization method of 

Figure 5. 

FIG. 8 is a perspective view of an optimized Na* • formaldehyde complex. 
FIG. 9 is a perspective view of an optimized formaldehyde dimer. 

Detailed Descrintlon nf thP invpntinn 
Embodiments of the presem invention will now be described with reference to the accompanying Figures, wherein like 
numerals refer to like elemems throughout. The tem^nology used in the description presented herein is not intended to be 
interpreted in any limited or restrictive manner, simply because it is being utilized in conjunction with a detailed description of 
certain specific embodiments of the invention. 

figure 1 is an illustration of a first group of atoms 10 and a second group of atoms 12 which are in a spatial 
relationship with one another. Atoms a and c of Rgure 1 are part of the first group 10 and atoms b and d are part of the second 
group 12. In some applications of the invention, each of the atoms in the first group 10 will be connected to at least one other 
atom of the first group 10 by a strong "bonded" interaction. Similarly, each of the atoms in the second group 12 will typically 
be connected to at least one other atom of the second group 10 by a strong bonded interaction. Interactions between atoms of 
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the first group 10 and the second group 12 are typically of a 'non-bonded" nature. Bonded interactions are characterized as 
strong, short range, and directional, while non-bonded interactions are generally characterized as weak, long range, and typically 
functions only of the distances r^ and r^ between non-bonded interacting pairs of atoms. 

It will be appreciated, however, that the atom groupings need not necessarily be based on the bonded vs. non-bonded 
distinction, and that the forces on any given atom may result from a combination of bonded and non-bonded interactions with 
any of the other atoms in either of the atom groups. In some appfications. the first group of atoms may comprise all of or a 
portion of a first molecule, and the second group of atoms may comprise aU or a portion of a second molecule. Alternatively, the 
first and second groups of atoms may comprise different portions of the same molecule. In other appDcations of the inventioa 
the second group of atoms comprises a reaction catalyst, and tiie first group of atoms is all or part of a substrate molecule. 
Interactions with surfaces, films, and crystal stmctures may also be defined as interactions between such first and second atom 
groups. 

The location and types of atoms present in the first and second groups of atoms define a potential field which 
determines the energy of a given configuration of the system. The potemial field may be expressed as a sum of temfis. with 
each l&m representing a particular class of interaction. Most force fields include temis which represent bonded interactions, 
which are also sometimes referred to as intramolecular Irteractions. Temis relating to bonded interactions may include 
configuration energy shifts from bond stretching, angle bending, out of plane defomatioa and possibly a variety of other 
molecular distortions from an unperturbed state. Any of the known parameterizations of the bonded interactions may be used 
with the present invertion. and these imramolecular potemial functions wiU not be described in furtiier detail herein. 

Non-bonded, or imemiolecular interaction potemials gcneraOy inchide a van der Weals interaction and a 
parameterization of electrostatic imeractions. The van der Waals temi may take a variety of fomis. and is atttactive at large 
atomic separations and repulsive at small atomic separations. Different known parameterizations of the van der Waals temi 
may be used with the present invention, and as with the intramolecular force parameters. Will not be discussed in greater detad 
herein. 

A general expression for a potential fieW used for molecular modeling may thus be described as foDows: 
(1 ) O - Intramolecular terms ■•• van der Waals tenns + electrostatic tenns 

The energy of a particular configuration of atoms is detemiined by summing the potential energy at points typically 
placed at atomic nuclear locations corresponding to the relevam atoms in the system being modeled. Forces acting on the 
atoms are determined by calculating the gradient of ti» potential energy field at the same atom centered pomts. 

Turning now to the electrostatic temis of the potential energy field, some models restrict tiieir consideration of 
electrostatic interactions to tiie assignment of a net charge to each atom. These net charges are used to compute a potential 
energy field at each atom centered point according to qiOj/r, for each atom pair, where q, is the partial charge of atom i, q, is tfie 
partial charge of atom j. and r, is the distimce between atom i and atom j. The DREIDING and CFF force fields described above 
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are examples of this type of field. It is one aspect of the invention that electrostatic interactions are treated in a more accurate 
manner. Specifically, each atom is assigned at least a dipole polarizabiiity tensor which may be anisotropic. In soma 
advantageous embodiments, each atom is modeled as a combination static charge, static dipole moment, static quadrupole 
moment, and a dipole polarizabiiity which, as mentioned above, may be anisotropic for at least some of the atoms of the 
modeled system 

In Figures 2 through 4, parameterization of the electrostatic interactions according to one embodiment of the invention 
is set forth. As is explained more fully below, the values of the above mentioned parameters are selected by ensuring that the 
overall electrostatic properties of molecules or other aggregates of parameterized atoms reproduce accurate quantum 
mechanically calculated electrostatic properties of those molecules or aggregates. In other words, the potential field and 
polarization response produced by an atom centered arrangement of the chosen charge, dipole moment, quadrupole moment, and 
polarizabiiity parameters will substantially match the potential field and polarization response produced by the actual motecuies 
or other aggregates of atoms that the model is intended to simulate. 

As will be appreciated by those in the art. the parameters associated with a given atom will be dependent on both the 
atom type and the other atoms it is bonded to. For example, an oxygen atom in an ether group will be parameterized differently 
15 than an oxygen atom in an aldehyde group. A large number of parameter sets may therefore be created which are devised to 

match those atomic aggregates which it is desired to model. If an interaction between a pair of molecules is to be modeled, as 
in an example set forth below, the atomic parameters are selected to reproduce molecular properties. Parameter derivations 
from selected functional groups rather than whole molecules may also be performed. 

Referring now to Figure 2, the creation of a parameter set for modeling the electrostatic interactions of a selected 
molecule begins at step 16. At this step, ab initio quantum mechanical calculations are used to compute values for a molecular 
dipole moment, molecular quadmpole momem, molecular octupole moment, molecular dipole polarizabiiity tensor, molecular 
quadrupole polarizabiiity tensor, and first derivatives with respect to atomic coordinates of ail of these molecular quantities. 
Producing these values for molecular multipole moments, polarizabilities, and their derivatives is computationally expensive, and 
is performed only once in the initial parameter set production process. 

At the next step 18. the desired atomic parameters are fitted to the previously determined molecular values. In one 
advantageous embodiment, the atomic parameters fitted to these quantum mechanical values include the atomic charges, 
dipoles, quadrupoles, and dipole polarizabilities discussed above. 

One set of atomic parameters which may be calculated at step 18 is the partial charge assigned to each atom. The 
assignmem of partial charges to atoms is a well established technique, and several methods of assignning charges to atoms in 
molecules or functional groups have been described. In one advantageous method, these values are arrived at by fitting tha 
charge parameters (and the derivative of the charges with respect to internal coordinates, also called a charge flux) to the 
calculated molecular dipole moment. The mathematical relationship of molecular dipole moments (and their derivatives) to 
atomic charge parameters in planar molecules has been previously described, and one suitable method to select atomic charge 
parameters so as to fit an ab initio molecular dipole moment is described in detail by U. Oinur and A.T. Hagler in J. Chem. Phys. 
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91. 2949-2958 11989). The disclosure of this J. Chem. Phys. article is incorporated herein by reference in its entirety. 

Static atomic dipole moment parameters are also fitted to the molecular characteristics. Molecular models 
incorporating static atomic dipole moments have also been produced. In one advantageous model, the static atomic dipole 
moment parameters are mathematically related to the molecular quadrupole moment and its derhfirtives. For planar molecules, 
this relationship, and the derivation of atomic dipole parameters using It is also described in the 1989 Dinur and Hagler article 
cited above. Step 18 may also comprise the calculation of static atomic quadnipoles. In one embodiment, these parameters 
may be fitted to the calculated molecular octupole moment and its derivatives. The constniction of permanent atomic 
quadrupole parameters from molecular octupole moments and their derivatives is described in U. Dinur, J. Comp. Chem. 12. 91- 
1 05 (1991). The disclosure of this J. Comp. Chem. article is incorporated herein by reference in its entirety. 

Additionally, step 18 may comprise a fit of atomic dipole polarizability tensor components to the molecular quadrupole 
polarizability tensor components and their derivatives. This procedure may be performed with manipulations which are 
described in additional detail below with reference to Figure 3. Following the derivation of the desired atomic parameter set at 
step 18. the method of Figure 2 moves to step 20 where interactions between groups of atoms are calculated using the atomic 

parameters in an atom-atom potential modd. 

With reference to the water molecule structure set forth in Figure 3. the deterniination of atomic dipole polarizability 
tensor components from molecular quadrupole polarizability tensor components and their derivatives will now be described. As a 
starting point use is made of the following expressions of the six unique elements of the molecular quadrupole moment in temis 
of atomic dipole moments and atomic partial charges: 



1 2 

2^' 'V' 



• L 

®» =Z [^"1.^/ -"ft-.V/ 
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Where 0,1, is the ab component of the molecular quadrupole moment. Uj, is the a component of the dipole moment of 
atom i, Qi is the partial charge of atom i, and x, y, and :„ are displacements in the x. y. and z directions of atom I From these 
equations, it is possible to express an induced molecular quadrupole moment in terms of induced atomic dipole moments as set 
forth below, if the explicit consideration of induced charges is ignored q- is set to 0): 

®;.=Z[2«,>,-«;jf,-«;r,] 

®>'=f Z i"iy^,+"'uy, ] 

In the above equations, induced moments are indicated with a prime. By considering Induced molecular quadrupole 
moments as sums of induced dipole moments caused by an external applied field, the molecular quadrupole polarizability tensor 
components can be related to components of atomic dipole polarizabiKty tensors for the atoms in the molecule. For example, in 
an applied field in the x direction, F,: 

I 

= ^xZ[2" -cc -a ,_^;;J 

I 

= ^x,a > therefore 

Where Ag, be is the a, be component of the molecular quardrupole polarizability tensor, and wherein ab is the ab 
component of the dipole polarizability tensor of atom i. Similar relationships may be derived in the same way by considering the 
other components of the induced quadrupole moment tensor and considering fields applied in the / and / directions. For 
example: 
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A.» =r [2a,..„Jc. -o-i^y, -^i^^i ] 

Z i 

A,^ =S [2a,,„*, -a,.^y, ] 

As these equations include sums over all the atoms in the molecule, there will be 6N unknowns, as there are six unique 
components of the dipole polarizabHity tensor for each of the N atoms in the molecule. Because there are only 18 unique and 
independent elements of the molecular quadrupole poiarizability tensor, there may be many more unknowns than independem 
equations. For many molecules and functional groups, therefore, this set of equations will not lead to a unique determination of 
the complete set of atomic dipole poiarizability parameters. For some classes of molecules, however, the additional 
consideration of the derhratiwes of the molecular quadrupole polarizabifity tensor components along with algebraic manipulation 
results in the simplification of the sums, and atomic dipole poiarizability parameters may be uniquely calculated from the 
derivatives of the components of the molecular quadrupole poiarizability tensor. 

As an example, for a planar molecule in the x-y plane, the derivative of the molecular quadrupole poiarizability tensor 
component Am, with respect to the z displacement of atom k is: 



This derwative pulls 0^, (the zz componem of the atomic dipole poiarizability tensor for atom k) out of the sum over 
all the atoms. As another example, the derivative of the molecular quadrupole poiarizability tensor componem A,^ with respect 
to the z dispbcement of atom k pulls o».„ (the xx component of the atomic dipole poiarizability tensor for atom k) out of the sum 
over all tiie atoms as seen below: 



In some cases, these equations can be combined to cancel out the remaining sums over all atoms. For Onear molecules 
located along the x-axis, where the y and z displacemem of each atom may be set to zero after perfomiing the derivatives, 
algebraic manipulation of the equations for the derivatives of the components of the molecular quadrupole poiarizability tensor 
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6L &4 dz^ 

oL dz^ 



For planar molecules, where only the z component of atomic 
derivatives, the following relationships may be derived: 



3[ &4 &4 



displacement can be set to zero after performing the 



a 



k.xy 



From the above equations, it can be seen that for linear molecules, the components of the atomic dipole polarlzability 
tensor for each atom can be computed arithmetically from the derivatives of the molecular quadnipole polarlzability tensor 
components, which are calculated, as mentioned above, with high accuracy ab initio quantum mechanical calculations. For 
planar molecules, the off-diagonal elements of the atomic dipole polarlzability tensor may be computed arithmetically, but only 
two equations exist for the three diagonal elements a,„, a,„, and For non-Bnear molecules therefore, the atomic dipole 
polarizability tensor components are not uniquely determined by the derivatives of the molecular quadnipole polarlzability tensor 
components alone. Accordingly, calculations for non-finear molecules require simplifications and/or infomiation derived from 
sources in addition to the molecular quadrupole polarlzability tensor derivatives. 

The computation of atomic dipole polarizability tensor components for the water molecule demonstrate this situation. 
In this case, coordinates can be selected such that the water molecule lies in the xy plane, with the oxygen atom at the origin as 
is illustrated in Figure 3. In this configuration: 
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Thus, o„, „ - o„i „ - J3xh. Using quantum mechanicaHy computed values for and x„. a value for 

ohi. a and om. a OSmi A' resuhs. 

The value for aa„ may be detemiined by additional consideration of the n component of the molecular dipoie 
polarizability for the vrater molecule because a»„« - ao^ + 2aHa. This leads to a value for ao.„ of 0.361 

The derivatives of the molecular quadrupole polarizability tensor of water may then be used acconling to the equations 
for planar molecules set forth above because the u component of the tensor has been specified, resulting in unique values for 
the other dipoie polarizability tensor components. When the quantum calculations of the molecular quadrupole polarizability 
tensor componem derivatives are detem«ned using HFI6.31 G* quantum mechanical calculations, the atomic dipoie polarizability 
tensor componems for the atoms in a vuatw molecule are (in units of A^: 

ao^ - 0.481, oo.„ - 0.561, - 0.361 



«H..» - ttHiB - 0.123, a„,,„ - Ohi„ - 0.221. a„,^ - a^M " 0-0367 
aHij>,-0-115,aH2.i,- -0.115 



As an isottopic dipoie polarizabifity requires an equal value for all diagonal elemems of the dipoie polarizability tensor 
and zeros for ail off diagonal elements, it may be noted that the atomic dipoie polarizabUity tensors for the oxygen and both 
hydrogens in the water molecule are anisotropic. Both the oxygen and hydrogens are more poiarizable in the y direction than 
they are in the x direction. Thus, for the water molecule as parameterized according to one embodiment of the imrentioa applied 
electric fields which are not exactly aligned with the x or y axis of the molecule wfll induce dipoie moments which are not 
parallel to the applied field. 

For molecules of N atoms m artjitrary three dimensional arrangements, the molecular quadrupole polarizability tensor 
components and their derivatives do not necessarily uniquely define the desired 6N atomic dipoie polarizability tensor 
components. It will therefore be appreciated that the use of molecular symmetries and infomation in addition to the molecular 
quadrupole moment tensor and its derivatives may be used to develop atomic dipoie polarizability tensors which fit molecular 
multipole quantities. In some cases, the atomic dipoie polarizabiDty tensor components can be fomwiated to fit the molecular 
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dipole polarizability tensor and its derivatives in addition to the molecular quardrupoie polarizability tensor and its derivatives. 
Furthermore, atomic parameters including the dipole polarizability tensor can be chosen to produce a simultaneous least squares 
fit to the multipole momems (and their derivatives) of a family of molecules, wherein each molecule of the family contains one or 
more atoms to be parameterized which are expected to behave similarly for each member of the family. It will also be 
appreciated that in many molecules and families of molecules, some of the temis in the equations above may be set to zero if 
their contribution is expected to be sufficiently smail. 

As illustrated in the flowchart of Figure 2. the general molecular modeling scheme is to parameterize a complex 
quantum mechanical system with atom centered force field parameters, and then to use those parameters to model the behavior 
of the complex system. A general modeling method using the parameters described above is presemed in Figure 4. The first 
step 28 of the method comprises parameterizing the behavior of a selected system with one or more of the parameters 
mentioned above. These parameters may include a charge, a static atomic dipole, a static atomic quadnipole, and an anisotropic 
atomic dipole polarizability. It will be appreciated that not all of the atoms in a system need be parameterized with the above 
described parameter set. Furthermore, one or more atoms of the system may be parameterized with only a subset of these 
parameters, and one or more of the atoms of the system may be parameterized with additional parameters not specifically 
15 described herein. 

In the next step 30. an electric field is calculated at the locations of one or more of the parameterized atoms of the 
system. This electric field may be due to the presence of the atoms, may be an applied field, or may be a combination of the 
two. In the presence of this electric field, induced dipole moments will be formed according to the atomic dipole polarizability 
parameters defined for the model. The assigned charges, static dipoles, induced dipoles. and static quadrupoles will also 
interact with the electric field local to each parameterized atom being considered. These effects are calculated at step 32. At 
step 34. one or more physical properties of the modeled system are predicted, based on the imeractions calculated at step 32. 
It wai be appreciated that a wide variety of physical properties may be determined in accordance with this embodiment of the 
invention. Optical properties, crystal structures, binding affinities, as well as other physical properties may be predicted. 
Calculated interaction energies may be used in several applications, including molecular dynamics simulations and Monte Cario 
calculations where the interaction energy at a variety of molecular conformations and positions is used. In addition, energy 
minimizations may beperfomied to determine stable geometries of molecular imeraction. 

hi the field of phamiaceutical research, the geometry and energy of ligand binding is of significant interest. This 
application of the invention is illustrated in Figure 5. Referring now to this Figure, one embodimert of a method for modeling 
geometries and energies of interaction between groups of atoms begins at a start state 40. It is assumed that the atom groups 
to be modeled are defined at this stage, and are placed in a starting configuration. In the phamiaceutical field, one group of 
atoms may comprise aV or a portion of a selected candidate molecule such as a ligand molecule, and one group of atoms may 
comprise a portion of a protein. In this application, it may be of interest to determine the binding energy and geometry of the 
ligand-protein complex. Aitemathrely, two atom groups may comprise different portions of a protein, and the modeling may be 
perfomied to evaluate an intramolecular non-bonded imeraction important in protein confomiation. Following a prediction of 
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physical properties using the mndellno techniques herein described, candidate molecules may be synthesized and tested for the 
predicted physical properties. 

It will be understood that in most applications, a general purpose computer is used to implement the methods 
described herein. The general purpose computer will have access to a data base of interaction parameters which is stored in a 
data storage device such as a CD-ROM. magnetic disk, semiconductor integrated circuit memory, or the like. The computer wiH 
also include a processor for accessing the data base, performing electric fieM calculatioms, computing induced dipole moments, 
calculating interaction energies, etc. Also inchided will be input and output devices for user interface such as keyboards, graphic 
display, printer, etc. The commands which configure such a general purpose computer to implement the methods of the 
imrentton are stored on a computer readable medium such as a CD-ROM for access by the general purpose computer. 

Returning to a discussion of the method of figure 5, at step 42, force field parameters appropriate to the atoms and 
atomic bonds of the system being modeled are retrieved from the data base. These parameters wiB typically include parameters 
related to bonded interactions, van der Weals parameters, and electrostatic interaction parameters for atoms in a wide variety 
of molecules and functional groups. Which parameters are retrieved will of course depend on the atoms and functional groups of 
atoms being modeled. 

At step 44, an electric field is calculated at the locations of the parameterized atoms. In tiiis embodiment, this electric 
field may be calculated using classical electrostatics based on the pointwise arrangement of atomic charges, static and induced 
dipoles. and quadrupoles defined by (1) the parameters retrieved at step 42 and (2) the coordinates of the parameterized atoms 
of the system. Because the fieW depends on the induced dipoles. and the induced dipoles depend on the field, the field and 
induced dipoles are n«de self consistent using a procedure described in more detaH with ref wence to Figure 6. 

Once the electric fieW is computed, at step 46 static multipole interaction energies and forces are determined. 
Classical electrostatics may be used to perform a palrwise computation of charge/charge, charge/static dipole. static 
dipote/static dipole, and charge/quadmpole interaction energies for the parameterized atoms of the system. At step 48. the 
polarization energy of the system is determined. The polarization energy is the contribution to the total energy which is 
attributable to the formation of mduced dipoles according to the retrieved atomic dipole polarizabilhy tensors. 

According to this embodiment of the invention, therefore, the electrostatic and polarization contribution of the 
interaction energy between the atoms of the system which is calculated at steps 46 and 48 of figure 5 may be expressed as: 



wherein: 



E^q - charge/charge interaction energy - ^QiQj I 
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Eqft, - charge/static dipole interaction energy - - • J^^, 



Eq/u- - charge/induced dipole interaction energy - - • /T* 

E„to - static dipole/static dipole interaction energy - (-l/2)2]w^ • JT" 

I 

E^u- - static dipole/induced dipole interaction energy - •-P" 
E„v - induced dipole/induced dipole interaction energy - - (l/2)^Uf • /T"' 
Eq/o - charge/quadrupole interaction energy = (J)^ 

E,rt - self polarization energy = (l/2)2]«; • (/Tj + + /T" ) 

'5 The total polarization energy contribution to the interaction energy is E^^,. + E„v E^- + E,.,,. which may be expressed 

as: 

(-i/2)Z[«,.ir;-f«;./r;] 
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Forces due to these interactions may be computed at steps 48 and 50 by perturbing the coordinates of the 
parameterized atoms and recalculating the interaction energies to determine the gradient of the interaction energy at the atomic 
coordinates. One process for determining the polarization energy gradient (and thus the polarization forces) is explained in more 
detail with reference to figure 7. At step 50, interaction energies and forces are computed for force field terms not described 
with reference to steps 46 and 48, including the mtramolecular and van der Waals energms and forces. 

If the current atomic configuration is a stable configuration of the system, the calculations performed at steps 46, 48, 
and 50 above will result in calculated forces on the atoms of near zero. At step 5^ therefore, it is determined whether or not 
these forces are below a threshold value. If they are not at step 54, the coordinates of the atoms of the system are 
incrememed. In one embodiment, the incrementing is doi» in the direction of the energy gradient although other energy 
minimization algorithms could be used with similar results. After the coordinates are incremented, the method loops back to 
step 44 to recalculate the field with the atoms at their new coordinates. An iterathre energy minimizatian process may therefore 
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be pert ormed until the forces on the atoms of the system are below the threshold at step 52. This indicates that a stable 
geometry has been reached, and at step 56 a computed total interaction energy and the stable atomic coordinates are output. 

Figure 6 is a flowchart illustrating an iterative method of making the local electric fields and the induced dipoles self 
consistent when the field is calculated at step 44 of Figure 5. The method begins at step 58, where the field due to the static 
5 charges assigned to each atom is detenrtned at each atomic location. Next at step 60. the field at each atomic location from 

the static dipoles and quadrupoles of each parameterized atom is calculated. Initially, at step 62, the induced atomic dipoles and 
the field produced by induced atomic dipoles for each atom are set to zero. At step 64, the induced dipoles formed by the 
presence of the current calculated field are determined using the anisotropic atomic dipole polarizability tensors discussed at 
length above. In the first iteration of the method of Figure 6, the field will be due solely to the contributions from the static 
1 0 multipoles calculated at steps 58 and 60, as the induced dipoles and the field resulting from them are initially zero. 

Next, at decision state 66. the new induced dipoles for each atom are compared to the old induced dipoles. If the 
difference for any of the parameterized atoms is greater than a threshold value, the dipoles and field are not considered self 
consistent. In this case, the field due to the induced dipoles is calculated at step 68, and the method loops back to step 64, 
where a new set of induced dipoles are determined. This new set of induced dipoles are computed with reference to the local 
1 5 fields as updated with the contribution from the induced dipoles calculated in the last iteration. 

This iterative loop continues until the new induced dipoles are sufficiently close to the last induced dipoles at step 66, 
which indicates self consistency between the electric field and the atomic dipole moments induced by the electric field. Once 
this point is reached, the field and induced dipoles are set to the values produced by the last field-induced dipole update iteration 
at step 70. Although an iterative method is illustrated in Figure 6. those in the art will recognize that an analytical matrix 
20 inversion method could also be used to accomplish the same result. 

Figure 7 illusuates an iteratwe method for calculating the polarization forces at step 48 of Figure 5. Referring back to 
the formula for the polarization energy given above, it can be seen that the gradient of the polarization energy may be written as: 

dE , „ a/r"' a/T^ sw' „g du. 

dXj i dxj dXj dXj dXj 

25 

Thus, a calculation of the polarization energy gradient requires the determination of the derivatives of the induced 
dipoles themselves, the derivatives of the field due to the induced dipoles, and the derivatives of the field due to the atomic 
charges. To determine the derivatives of the induced dipoles. it is additionally required to determine the derivatives of the field 
due to the static atomic dipoles. In an iterative manner analogous to that illustrated in Figure 6, the derivatives of the field due 
30 to the induced dipoles is made self consistent with the derivatives of the induced dipoles themselves. This method begins at 

step 74 of Figure 7 where the derivatives of the field due to the atomic charges and the field due to static dipoles are 
determined. Next, at step 76. the derivatroes of the induced dipoles and the field due to the induced dipoles is set to zero. At 
step 78, the derivatives of the induced dipoles are then calculated. In the first iteratioa tiiis derivative is calculated using only 
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the derivatives of the fields due to the charges and the static dipoles, as the derivatives of the induced dipoles is set to zero. 

The method next moves to a decision state 80, where the new values of the derivatives of the induced dipole 
moments are compared to the prior derivatives of the induced dipole moments for each parameterized atom. If the result of 
these comparisons is more than a given threshold, the induced dipole derivatives are not consistent with the induced dipole field 
5 derivatives. In this case, the method moves to step 82 and new derivatwes of the induced dipole field are calculated using the 

derivatives of the induced dipoles calculated in the last iteration at step 78. From step 82, the method loops back to step 78, 
and recalculates new derivatives of the induced dipole moments using the derivatives of the induced dipole field derivatives 
calculated in step 82. 

If at step 80 it is determined that the differences between the newly calculated derivatives of the induced dipole 
10 moments and the prior derivatives of the induced dipole moments from the last iteration are less than a threshold value, the 

induced dipole derivatives are considered consistent with the induced dipole field derivatives. In this instance, the method moves 
to step 84, where the polarization energy gradient is determined using the final values for the induced dipole derivatives and 
induced dipole field derivatives. 

With reference now to Figures 8 and 9, examples of interaction energies and geometries calculated using a force field 
15 including static atomic multipoles and anisotropic atomic dipole polarizabilities are presented. In the first example, illustrated in 

Figure 8, the structure and interaction energy of the formaldehyde/sodium ion complex was modeled. The force field used in this 
model included charge, static dipole, and static quadruple parameters, as well as atomic dipole polarizability tensor components 
as described above. The van der Waals parameters were adjusted such that the optimized Na*/0 distance 90 equaled 2.16 A, 
which is the value calculated using ab initio quantum mechanical methods. 
2" The interaction energies determined by the model (in kcal/mole) are: 



intramolecular: 0.02 

van der Waals: 7.88 

charge/charge: -25.40 

25 charge/static dipole: .2.61 

charge/static quadrupole: 2.77 

polarization: -12.29 

TOTAL: -29.83 



These values may be compared to a quantum mechanically calculated total interaction energy of -28.3 kcal/mole 
(generated using HF/6-31G* with a basis set superposition enror correction). In contrast, using a force field model which 
includes only charge/charge electrostatic interaction parameterization, the polarization energy is neglected, and a value of 
approximately -18 kcal/mole is obtained. The inclusion of atomic dipole polarizabilities in the model parameter set thus greatly 
improves the match between the predictions of the force field calculations and the computationally expensive quantum 
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nriechanical calculations. 

Referring now to Figure 9. an optintized structure of the fomraldehyde dimer is illustrated. Using quantum mechanical 
calculations |MP2;6-31 +G(2df, 2pd)), the predicted optimized structure comprises a perpendicular stacked structure, with the 
distance between atoms 02 and H7 (denoted 92 in Figure 9) being 2.38 A. and the distance between atoms 06 and H3 (denoted 

5 94 in Figure 9) being 2.81 A. The quantum mechanically calculated interaction en»gy is 4.6 kcal/mole. 

Using a force field including charge, static dipole, static quadrupole parameters, and atomic dipole polarizabiBty tensor 
components, the optimized structure produced by the model is the perpendicular stacked structure wherein the distance 
between atoms 02 and H7 92 is 2.77 A, and the distance between atoms 06 and H3 94 is 3.05 A. The calculated interaction 
energy using this parameterized model is 4.5 kcal/mole. 

1 0 If the parameterized model includes only charge/charge electrostatic interaction terms, the wrong optimized structure 

is predicted. The upper formaldehyde molecule lies flat on top of the lower formaldehyde molecule, and the 02 H7 and 06 H3 
distances are both equal to 3.17 A. By running the formaldehyde dimer optimization with differem models which consider only 
subsets of the charge/charge, charge/static dipole, static dipole/static dipole, charge/quadrupole, and polarization interactions, it 
has been found that the charge/quadrupole interaction appears significant in the formation of the optimized perpendicular 

1 5 stacked structure rather than the flat stacked structure. 

Molecular models in accordance with the invention therefore provide increased accuracy in predicting the physical 
properties of subject molecular systems. This will be especially true in those systems which include significant anisotropic 
polarizations, such as systems including ion-benzene interactions. Because the parameter set is created to fit higher order 
molecular moments and polarizabilities, more accurate energy calculations should result, as these molecular properties often 

20 contribute significantly to interaction energies, and are not adequately taken into account in current modeling schemes. 

Furthermore, the ratio of parameters to quantum observables can be increased by fitting, for an N atom molecule, 6N atomic 
dipole polarizability parameters to 18 molecular quadmpole polarizability tensor components and 54N derivatives of those 
components. 

The foregoing description details certain specific embodimems of the invention. It will be appreciated, however, that 
25 no matter how detailed the foregoing appears in text, the invention can be practiced in many ways. As is also stated above, it 

should be noted that the use of particular temrinology when describing certain features or aspects of the invention should not be 
taken to imply that the tenninology is being re-defined herein to be restricted to including any specific characteristics of the 
features or aspects of the invention with which that terminology is associated. The scope of the invention should therefore be 
constnied in accordance with the appended Claims and any equivalents thereof. 

30 
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WHAT IS CLAIMED \Rr 

1. A computer implemented method of modeling an interaction between a first group of one or more atoms 
and a second group of one or more atoms, wherein at least a portion of said first group of one or more atoms is 
characterized by a quadrupoie polarizability tensor, said method comprising: 
5 defining, for at least a first one of said first group of atoms, elements of an atomic dipole polarizability 

tensor utilizing derivatives of elements of said quadrupoie polarizability tensor; 

calculating a dipole moment induced in said first atom from a local electric field, said local electric field 
being dependent on a relative orientation and separation of said first group of one or more atoms and said second 
group of one or more atoms; 

determining an interaction energy between said first group of one or more atoms and said second group 
of one or more atoms which includes a contribution from said induced dipole moment; and 

outputting said interaction energy. 
2. The method of Claim 1, wherein said first group of one or more atoms comprises a first molecule, and 
said second group of one or more atoms comprises a second molecule. 
^5 3. The method of Claim 1, wherein said first group of one or more atoms comprises a first portion of a 

molecule, and wherein said second group of one or more atoms comprises a second portion of said molecule. 

4. The method of Claim 1, wherein said second group of one or more atoms comprises a catalyst 
5» The method of Claim 1, wherein said atomic dipole polarizability is anisotropic. 
6. The method of Claim 5, wherein said electric field and said induced dipole moment are not parallel. 
20 7. The method of Claim 1, wherein said first molecule is also characterized by a molecular quadrupoie 

moment tensor, and wherein said method additionally comprises defining elements of an atomic dipole moment tensor for 
said first atom utilizing derivatives of said molecular quadrupoie moment tensor. 

8. The method of Claim 1, wherein said first molecule is also characterized by a molecular octupole moment 
tensor, and wherein said method additionally comprises defining elements of an atomic quadrupoie moment tensor for said 

25 first atom utilizing derivatives of said molecular octupole moment tensor. 

9. The method of Claim 8, wherein said calculated interaction energy additionally comprises a contribution 
from said atomic quadrupoie. 

1 0- A method of evaluating a molecule for suitability for a particular purpose comprising: 
selecting a candidate molecule comprising a plurality of atoms; 
^0 calculating a dipole moment induced in at least one of said plurality of atoms from a local electric field, 

said induced dipole moment being non-parallel to said local electric field; 

predicting one or more physical properties of said candidate molecule using said calculated dipole 

moment. 
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11. The method of Claim 10, wherein predicting one or more physical properties of said candidate molecule 
comprises calculating an interaction energy between said candidate molecule and a target molecule which includes a 
contribution from said induced dipole moment. 

12. The method of Claim 11, wherein said interaction energy additionally includes a contribution from a 
permanent quadrupole moment in said first atom. 

13. The method of Claim 12, wherein said interaction energy additionally includes a contribution from a 
permanent dipole moment in said first atom. 

14. The method of Claim 13, wherein said interaction energy additionally includes a contribution from a 
charge on said first atom. 

1 5. The method of Claim 1 0, additionally comprising: 
synthesizing said candidate molecule; and 

testing said synthesized candidate molecule for said one or more physical properties. 

16. The method of Claim 10, wherein said induced dipole moment calculation comprises using an atomic 
polarizability tensor which may be anisotropic. 

17. The method of Claim 16, wherein said candidate molecule is characterized by a molecular quadrupole 
polarizability tensor, and wherein the elements of said anisotropic atomic dipole polarizability tensor have been determined 
utilizing derivatives of elements of a molecular quadrupole polarizability tensor. 

18. A method of evaluating a candidate molecule for suitability for a particular purpose, said candidate 
molecule comprising a plurality of bonded atoms, said method comprising: 

parameterizing behavior of said candidate molecule with a plurality of atomic parameters associated 
with at least one of said plurality of bonded atoms, said plurality of atomic parameters including elements of an 
anisotropic atomic dipole polarizability tensor; 

determining a dipole moment induced in said at least one of said plurality of atoms due to a local electric 
field using said atomic dipole polarizability tensor; and 

predicting one or more physical properties of said candidate molecule using said induced dipole moment. 

1 9. The method of Claim 1 8, additionally comprising: 
synthesizing said candidate molecule; and 

testing said synthesized candidate molecule for said one or more physical properties. 

20. The method of Claim 18, wherein predicting one or more physical properties of said candidate molecule 
comprises calculating a polarization energy. 

21. The method of Claim 18, wherein predicting one or more physical properties of said candidate molecule 
comprises determining an energy of interaction between said candidate molecule and a target molecule. 

22. The method of Claim 18, wherein said plurality of parameters additionally includes a permanent atomic 
quadrupole moment. 
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23. The method of Claim 22, wherein said plurality of parameters additionally includes a permanent atomic 
dipole moment. 

24. The method of Claim 23, wherein said plurality of parameters additionally includes an atomic charge. 

25. The method of Claim 24, additionally comprising calculating an energy of interaction between said local 
electric field and said permanent atomic quadrupole moment, permanent atomic dipole moment, and atomic charge, and 
wherein predicting one or more physical properties of said candidate molecule comprises using said energy of interaction 
between said candidate molecule and said target molecule. 

26. A computer readable medium having stored thereon commands which cause a general purpose computer 
to perform a method of modeling interactions between a first group of one or more atoms and a second group of one or more 
atoms, said method comprising: 

retrieving, for at least one of said first group of atoms, elements of an anisotropic atomic dipole 
polarizability tensor from a data storage device; 

modeling an electric field which is dependent on a relative orientation and separation of said first group 
of atoms and said second group of atoms; 

calculating a dipole moment induced in said at least one of said first group of atoms from said electric 
field using said anisotropic polarizability tensor; 

calculating an interaction energy between said first group of atoms and said second group of atoms 
which includes a contribution from said induced dipole moment; and 

outputting said interaction energy. 

27. The computer readable medium of Claim 26, wherein said method additionally comprises retrieving, for 
said at least one of said first group of atoms, a permanent atomic quadrupole moment, a permanent atomic dipole moment, 
and an atomic charge. 

28. The computer readable medium of Claim 27, wherein said method additionally comprises calculating an 
interaction energy between said first group of atoms and said second group of atoms which includes a contribution from 
said permanent atomic quadrupole moment, permanent atomic dipole moment, and atomic charge. 

29. An apparatus for modeling the geometry and energy of interaction between first and second groups of 
atoms comprising: 

a data storage device storing an anisotropic atomic dipole polarizability tensor for at least one of said 
first group of atoms; 

a processor which (1) models an electric field produced at least in part by said first and second groups of 
atoms, (2) retrieves said anisotropic dipole polarizability tensor, (3) calculates a dipole moment induced in said at 
least one of said first group of atoms by said electric field, and (4) calculates an interacrion energy between said 
first and second groups of atoms which includes a contribution from said induced dipole moment; and * 
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an output device which reports said calculated interaction energy or other physical property derived 
therefrom. 

30. The method of Claim 29, vuherein said memory additionally stores a static atomic quadruple moment, a 
static atomic dipole moment, and a charge for said at least one of said first group of atoms, and w/herein said interaction 

5 energy includes a contribution from said static atomic quadrupole moment, said static atomic dipole moment and said 

charge. 

31. A method of modeling a stable geometry of interaction between first and second groups of atoms 
comprising: 

parameterizing electrostatic behavior of said first and second groups of atoms with a plurality of atomic 
10 parameters associated with at least one of said first group of atoms, said plurality of atomic parameters including 

an anisotropic atomic dipole polarizability tensor; 

calculating dipole moments induced in said at least one of said first group of atoms using said atomic 
dipole polarizability tensor; 

calculating an interaction energy between said first and second groups of atoms, said interaction energy 
1 5 including a contribution from said induced dipole moment; and 

iteratively optimizing a relative orientation of said first and second groups of atoms so as to reduce said 
calculated interaction energy. 

32. The method of Claim 31, wherein said first group of one or more atoms comprises a first molecule, and 
said second group of one or more atoms comprises a second molecule. 

20 33. The method of Claim 31, wherein said first group of one or more atoms comprises a first portion of a 

molecule, and wherein said second group of one or more atoms comprises a second portion of said molecule. 

34. The method of Claim 31, wherein said plurality of parameters additionally comprises a charge, a 
permanent dipole moment, and a permanent quadrupole moment for each of said plurality of atoms. 

25 
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